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■ We use very large cosmological N-body simulations to obtain accurate predictions 

C^) for the two-point correlations and power spectra of mass-limited samples of galaxy 

clusters. We consider two currently popular cold dark matter (CDM) cosmogonies, a 
critical density model (tCDM) and a flat low density model with a cosmological con- 
stant (ACDM). Our simulations each use 10 9 particles to follow the mass distribution 
O within cubes of side 2/i _1 Gpc (rCDM) and 3/i _1 Gpc (ACDM) with a force resolution 

I ' better than 10 -4 of the cube side. We investigate how the predicted cluster correla- 

tions increase for samples of increasing mass and decreasing abundance. Very similar 
behaviour is found in the two cases. The correlation length increases from rrj = 12 
13/i _1 Mpc for samples with mean separation d c = 30/i _1 Mpc to ro = 22 - 27/i -1 Mpc 
for samples with d c — 100ft. - x Mpc. The lower value here corresponds to rCDM and the 
upper to ACDM. The power spectra of these cluster samples are accurately parallel to 
those of the mass over more than a decade in scale. Both correlation lengths and power 
spectrum biases can be predicted to better than 10% using the simple model of Sheth, 
C$ ' Mo & Tormen (2000). This prediction requires only the linear mass power spectrum 

and has no adjustable parameters. We compare our predictions with published results 
for the APM cluster sample. The observed variation of correlation length with richness 
agrees well with the models, particularly for ACDM. The observed power spectrum 
(for a cluster sample of mean separation d c = 31ft~ 1 Mpc) lies significantly above the 
predictions of both models. 

Key words: cosmology:theory - dark matter - gravitation - galaxy clusters- 
simulations 



1 INTRODUCTION 

The last two decades have established cosmological N-body 
simulations as the principal tool for studying the evolution 
of large-scale structure. The earliest systematic studies used 
10 3 to 2 x 10 4 particles to follow evolution from white noise 
or other similarly ad hoc initial conditions (Gott, Turner & 
Aarseth 1979; Efstathiou & Eastwood 1981). They showed 
that nonlinear growth could produce a power law autocor- 



relation function similar to that measured for galaxies. Soon 
thereafter, the suggestion that the dark matter might be a 
weakly interacting massive particle led to the first simula- 
tions from initial conditions based on a detailed treatment 
of the physics of earlier evolution. These represented the 
dark matter distribution within cubic regions with periodic 
boundary conditions using only 3 x 10 4 particles. They were 
nevertheless able to show that, for adiabatic fluctuations 
produced during inflation, a neutrino-dominated universe is 
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not viable (White, Frenk & Davis 1983) while a cold dark 
matter (CDM) dominated universe is much more promising 
(Davis et al. 1985). 

Since this early work, many groups have used their own 
simulations to compare predictions for large-scale structure 
with the wealth of data coming from new observational sur- 
veys. As algorithms and computers have improved, the num- 
ber of particles treated in high resolution simulations has 
increased. Thus, White et al. (1987a, b) could already use 
2.6 x 10 5 particles to study CDM universes, while Warren et 
al. (1992), Gelb & Bertschinger (1994), Jenkins et al. (1998) 
and Governato et al. (1999) studied large-scale structure us- 
ing high resolution simulations with Ixl0 6 ,3xl0 6 , 1.7xl0 7 
and 4.7x 10 7 particles, respectively. More particles are better 
for two reasons. One can choose to have better mass reso- 
lution so that the internal properties of each structure are 
better defined, and one can simulate larger volumes so that 
more structures are included and the statistical distribution 
of their properties is better defined. Here, we report results 
for the spatial distribution of galaxy clusters from simula- 
tions using 1 x 10 9 particles. The volumes simulated are 
much larger than any attempted previously and are large 
compared even to the biggest currently planned observa- 
tional surveys. As a result our theoretical predictions have 
high precision. 

The two-point correlation function of galaxy clusters 
has been controversial for decades. The early work of Hauser 
& Peebles (1973) showed that rich galaxy clusters are more 
strongly clustered than galaxies, and estimates of the auto- 
correlation function of Abell clusters by Bahcall & Soneira 
(1983) and Klypin & Kopylov (1983) agreed on a power-law 
form which parallels the galaxy autocorrelation function but 
with substantially greater amplitude. Subsequent work has 
failed to agree on the strength of this enhancement and on 
its dependence on the properties which define the cluster 
sample. Thus, Sutherland (1988) argued that much of the 
apparent clustering in the original samples was induced arti- 
ficially by Abell's criteria for defining clusters. This conclu- 
sion has been supported by some subsequent studies (e.g. 
Croft et al. 1997 and references therein) and disputed by 
others (e.g. Olivier et al. 1993 and references therein). 

Although all authors agree that richer clusters are more 
strongly clustered, the strength of this trend is also disputed. 
Bahcall and co-workers (e.g. Bahcall & Cen 1992, Bahcall 
& West 1992) have argued that the correlation length, ro, 
defined via £ c ;(ro) = 1, scales linearly with intercluster sep- 
aration, rf c , 

r = 0.4d c = 0.4n c - 1/3 , (1) 

where n c is the number density of clusters above the cho- 
sen richness threshold. This scaling might be expected in 
a fractal model of large-scale structure (Szalay & Schramm 
1985) and appeared consistent with early measurements for 
Abell clusters (e.g. Bahcall & Soneira 1983). Other work has 
suggested that this apparent scaling reflects incompleteness 
in the Abell samples (e.g. Efstathiou et al. 1992, Peacock 
& West 1992). The more objectively defined APM cluster 
sample appears to show a significantly weaker trend of clus- 
tering strength with richness (Efstathiou 1996; Croft et al. 
1997). 

Quite surprisingly, both camps have used N-body sim- 
ulations of standard CDM cosmogonies to support their 



views. Bahcall & Cen (1992) found ro to increase roughly 
in proportion to d c for their simulated clusters, while Croft 
& Efstathiou (1994) found a weaker dependence. The latter 
authors found their cluster correlation function to be insen- 
sitive to the cosmic matter density, fl, and to depend weakly 
on the normalization of the power spectrum, as, but strongly 
on its shape. (Here, as denotes linearly extrapolated present- 
day rms mass fluctuation in spherical top hat spheres of ra- 
dius Sh' 1 Mpc .) Similar conclusions were reached by Eke et 
al(1996) who studied systematics in simulated cluster cor- 
relations, in particular the influence of the definition of a 
cluster. They argued that the different scalings of ro with 
d c seen in previous N-body simulations stemmed primarily 
from the use of different algorithms to identify clusters in the 
simulations. All this work suffered from the relatively small 
volumes simulated, which limited the statistical accuracy 
of the correlation estimates, especially for rare and massive 
clusters. A substantial improvement came with the work of 
Governato et al. (1999) who used more particles and treated 
significantly larger volumes. The simulations presented be- 
low provide a further major improvement by using 20 times 
as many particles and increasing the volumes treated by 
about two orders of magnitude. 

The results we present below are in general agree- 
ment with those of Governato et al. (1999), but our work 
achieves substantially higher statistical precision. We show 
very clearly that the strength of cluster correlations is pre- 
dicted to increase significantly with cluster richness in cur- 
rently popular CDM cosmogonies. Furthermore, these cor- 
relations can be predicted remarkably accurately (and with 
no free parameters) by the recent analytic model of Sheth, 
Mo & Tormen (1996). In this model, which refines that of 
Mo & White (1996), the two-point correlation function of 
dark haloes is proportional to that of the dark matter, the 
ratio of the two depending on halo mass and on the linear 
power spectrum of mass density fluctuations (see below). Mo 
& White tested their original model on a set of scale-free N- 
body simulations, finding good qualitative agreement. For 
CDM models, Sheth, Mo & Tormen found the quantitative 
prediction both of halo mass functions and of halo correla- 
tions to be improved substantially by generalising the Mo 
& White approach to ellipsoidal (rather than spherical) col- 
lapse. Our results here reach higher precision and extend 
these tests to rarer objects; a preliminary account was pub- 
lished in Colberg et al. (1998), which is superseded by the 
current paper. 

The second order statistics of the spatial distribution 
of clusters can, of course, be analysed using power spectra 
rather than correlation functions. Such an approach is par- 
ticularly advantageous for analysing fluctuations on large 
spatial scales. Recent observational estimates of the cluster 
power spectrum have been given by Borgani et al. (1997) and 
Retzlaff et al. (1998), and by Tadros, Efstathiou & Dalton 
(1998) for the Abell- ACO and APM clusters, respectively. 
For both samples, there is an indication of a peak in the 
power spectrum at a wavenumber of k ~ 0.03 h Mpc -1 . This 
is roughly coincident with the scale where a peak is expected 
for currently popular CDM models. The simulation data we 
present below verify that the cluster power spectrum should 
indeed parallel that of the mass on these scales. 

In the following Section we briefly discuss the Hubble 
Volume simulations and the way we have defined cluster 
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Table 1. Parameters of the Hubble Volume simulations 

Model H A h T cr 8 L bo x 
rCDM 1.0 0.0 0.5 0.21 0.6 2000 Mpc/h 
ACDM 0.3 0.7 0.7 0.17 0.9 3000 Mpc/h 




samples within them. In Section we present two-point 
correlations for these samples and compare them with the 
analytic model. In Section ^, we present power spectra for 
samples constructed to correspond directly to the APM clus- 
ter survey; an interesting result is that the observations and 
predictions are significantly discrepant for the current "best 
buy" cosmogony. We conclude with a summary of our main 
results. 



2 CLUSTERS IN THE HUBBLE VOLUME 
SIMULATIONS 

The two simulations analysed in this paper were carried out 
in 1997 and 1998 on 512 processors of the CRAY T3E at the 



Garching Computer Centre of the Max Planck Society. They 
used a specially stripped down version of parallel Hydra, the 
workhorse code of the Virgo Supercomputing Consortium. 
Details may be found in MacFarland et al. (1998). This code 
maximises the efficiency of memory use on the machine and 
allowed the trajectories of 10 9 particles to be integrated ac- 
curately, with a gravitational force resolution of about 10 -4 
of the side of the computational volume. Each simulation 
used about 50,000 processor hours of CPU time. The two 
cases studied were a (2000/i _1 Mpc) 3 volume of a rCDM uni- 
verse and a (3000/!. _1 Mpc) 3 volume of a ACDM universe. In 
both cases the mass of a single particle is 2 x 10 12 hT 1 Mq 
and the simulation is normalized to yield the observed abun- 
dance of rich clusters at z — (White, Efstathiou & Frenk 
1993; Eke, Cole & Frenk 1996). These normalisations are 
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Figure 2. The two-point correlation functions of the tCDM (lower plots) and ACDM (upper plots) models for d c = 50fc _1 Mpc. This 
figure compares results from the simulations (dots with errorbars) with the linear (dashed line) and nonlinear (solid line) predictions 
from eqn M with the SMT prediction for b. For the ACDM model all quantities have been shifted upwards by one order of magnitude. 
la errorbars are plotted, as in fig. [f| 



also consistent with the level of fluctuations measured by 
COBE. The parameters of the simulations are summarized 
in table |IJ (Here, F denotes the spectral shape parameter; 
c.f. Efstathiou, Bond & White 1992.) 

The Hubble Volume simulations are essentially larger 
realisations of two of the cosmological models previously 
simulated by Jenkins et al (1998). We have checked that 
basic properties of our new simulations, such as the mass 
power spectrum and the velocity field, are consistent with 
expectations based on our smaller simulations. High order 
clustering statistics in the rCDM Hubble Volume simula- 
tion have been extensively studied by Szapudi et al (2000) 
and Colombi et al (2000) . Both Hubble Volume simulations 
were used by Jenkins et al (2000) in a study of the mass 
function of dark matter halos. The mass functions from the 
Hubble simulations are consistent with those from smaller 
simulations in the regions of overlap. 

Clusters of galaxies were identified in these simulations 
using a spherical overdensity (SO) group finder (Lacey & 



Cole 1994). This defines the cluster boundary as the sphere 
within which the mean density is 180 and 324 times the 
critical value in the rCDM and ACDM cases respectively. 
The lowest mass clusters considered in our analysis have 
75 and 39 particles respectively in the tCDM and ACDM 
models. We have checked that our results in the form we 
present below are insensitive to this choice. For example, we 
obtain almost identical results if clusters are defined using 
a friends-of-friends algorithm (Davis et al. 1985) with link- 
ing lengths of 0.2 and 0.164 in each model (which produces 
clusters with at least 86 and 44 particles respectively.) The 
choice of grouping algorithm and associated parameters af- 
fects the masses assigned to clusters in a systematic way, 
but has no significant systematic effect on their positions or 
on their ranking in mass. 

We construct a series of mass-limited cluster catalogues 
and characterise each one by the mean separation d c of the 
clusters it contains. The advantage of this parameterisation 
is that it allows a precise comparison with observed richness- 
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Figure 3. Correlation length, ro, as a function of mean intercluster separation, d c , for the tCDM (open squares) and ACDM (filled 
squares) simulations. The predictions of the SMT model are shown as solid lines. Also shown are data from the APM cluster catalogue 
(open triangles), taken from Croft et al. (1997). 



limited samples and with analytic models without any need 
to ensure that the mass definitions in the three cases cor- 
respond exactly. In the current paper, we consider only the 
clustering of clusters at z = 0. The excellent fits we find to 
the analytic theory also hold at other redshifts, so in practice 
the analytic formulae can be used to describe supercluster- 
ing at any redshift. 



3 TWO-POINT CORRELATION FUNCTIONS 
3.1 The Analytic Model 

Starting from a "Press-Schechter" (1974) argument similar 
to those in Bond et al. (1991) and Lacey & Cole (1993), Mo 
& White (1996) developed an analytic theory for the spatial 
clustering of dark haloes in hierarchical clustering models 
such as the many variants of CDM. They find that the two- 
point correlation function of dark matter haloes of mass M 
may be approximately related to that of the mass by 



&(r;Af) = b 2 (M)£(r 
where 
b{M) = 1 + 



<7 2 (M) 



(2) 



(3) 



Here, S c = 1.686 is the interpolated linear overdensity at 
collapse of a spherical perturbation, and a(M) is the rms 
linear fluctuation in overdensity within a sphere which on 
average contains mass M. Notice that although a(M) can 
be calculated directly from linear theory, £(r) in eqn |^ is 
the full nonlinear correlation function of the mass density 
field. This can be estimated from the linear-theory power 
spectrum using, for example, the approximation of Peacock 
& Dodds (1996). Thus, the nonlinear correlation function of 
haloes can be predicted without the need to carry out an N- 
body simulation. As shown by Cole & Kaiser (1989), eqn ^ 
can be derived by calculating how the Press-Schechter mass 
function responds to small changes in the threshold 8 C . 

It has long been known that the Press-Schechter mass 
function is not a perfect match to the mass functions found 
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Figure 4. The upper panel shows power spectra for galaxy clusters with d c = 30.9 h _1 Mpc from the rCDM simulation (filled dots), for 
the APM cluster sample (triangles; taken from Tadros et al. 1998), and for the dark matter in the simulation (open squares). The lower 
panel gives a bias factor defined as the square root of the ratio of the cluster and dark matter power spectra. The horizontal dotted line 
is the value of this bias predicted by the SMT model. 



in simulations (e.g. Efstathiou et al. 1988), and recent work 
has demonstrated that there is a corresponding system- 
atic shift in the bias calculated using the Cole-Kaiser argu- 
ment (Jing 1998; Sheth & Tormen 1999). Sheth, Mo & Tor- 
men (2000; SMT) have shown how the inclusion of a mass- 
dependent absorbing barrier in the excursion set deriva- 
tion of the mass function (Bond et al. 1991) can model 
the anisotropic collapse of cosmic structure and substan- 
tially improve the agreement between analytic theory and 
numerical simulation. Following the logic of Mo & White's 
extension of the excursion set formalism but using this mass- 
dependent threshold, SMT predict halo clustering in good 
agreement with simulation data. For our purposes, the effect 
of the SMT revision is to predict a slightly different b(M) 
from that in eqn. (3). 

A technical problem arises when comparing such ana- 
lytic formulae with simulations; it is unclear how to define 
the boundaries of simulated clusters so that their mass cor- 
responds to the mass M in eqn ^[ Although this might seem 



to introduce an additional degree of freedom, we can elim- 
inate it by using the corresponding analytic expression for 
the abundance of clusters to convert from sample limiting 
mass, M, to mean cluster separation d c . The predicted cor- 
relations can then be compared to those of a mass-limited 
sample of simulated clusters with the same mean separa- 
tion. This comparison has no adjustable parameters. Note 
that for such mass-limited samples the factor b 2 in eqn ^ 
is the square of the mean bias obtained by weighting b(M) 
by the abundance of clusters of that mass (see, for example, 
Baugh et al. 1998 and Governato et al. 1999). 



3.2 Results 

Figure Q shows cluster correlation functions for the rCDM 
simulation for mass-limited samples of clusters with mean 
separations of 40, 70, 100 and 130 h~ 1 Mpc. These samples 
contain 125,000, 23,000, 8,000 and 3,600 clusters respec- 
tively. We have computed la errors from the numbers of 
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Figure 5. As in fig. [l] but for the ACDM model. For this model the SMT prediction of the bias factor is high but only by 6%. 



pairs in each separation bin. Clearly, more massive clusters 
are more strongly clustered. Note also the very small error 
bars on these correlation estimates which are a consequence 
of the very large volume of our simulations. 

Figure H shows the correlation functions of samples with 
dc = 50/i~ Mpc from our two simulations, together with 
predictions from the SMT model of the last subsection. The 
predictions are shown separately for the two cases where 
£(r) is simply taken as the Fourier transform of the linear 
power spectrum, and where it is estimated using the non- 
linear model of Peacock & Dodds (1996). The correlation 
functions are very similar in the two cosmologies, showing 
that the strength of superclustering is not a good estima- 
tor of for CDM models normalised to match the observed 
abundance of clusters and having a mass correlation func- 
tion with a similar shape to the galaxy correlation function 
on large scales. (Note that the curves for ACDM have been 
raised by an order of magnitude for clarity.) The analytic 
predictions are in excellent agreement with the numerical 
results, particularly for £a ~ 1. Over the relevant range 
of scales the linear and nonlinear predictions for £(r) are 



quite close, and using the nonlinear formula gives at best a 
marginal improvement in the fit to the simulation results. 

In figure ^ we quantify the increase in clustering 
strength with cluster mass by plotting the correlation length, 
ro, of our mass-limited cluster samples as a function of their 
mean intercluster separation, d c . We estimate correlation 
lengths from plots similar to those of figure [l] by interpolat- 
ing between the points on either side of £ fe = 1. This figure 
again shows that our simulated volumes are large enough 
to estimate correlation lengths with high accuracy. The val- 
ues of ro for ACDM exceed those for rCDM by between 
10 and 20%. In both models, the increase in ro with d c is 
quite strong, although weaker than the direct proportional- 
ity suggested by Szalay & Schramm (1985) and Bahcall & 
Cen (1992). For rCDM the analytic prediction of ro is ac- 
curate to within a few percent on all scales; for ACDM it is 
about 10% high. 

The same general trend of ro with d c is also apparent in 
the simulations of Governato et al (1999) who considered an 
fi = 0.3 open CDM model (OCDM) and an tt = 1 standard 
CDM model (SCDM) . The clustering amplitude of clusters 
in OCDM is similar to that in ACDM, while that in SCDM, 
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although qualitatively similar, has much lower amplitude, 
reflecting the relatively small amount of large-scale power 
in this model compared to the other three. 

Our predictions may be readily compared with the 
measured values of ro for APM clusters given by Croft et 
al. (1997). Comparison with these data is relatively sim- 
ple because this cluster sample is approximately volume- 
limited. By contrast, comparison with X-ray selected clus- 
ter samples (e.g. Ebeling et al 1996, Guzzo et al 1999), 
which are flux-limited, requires more extensive modelling 
(see Moscardini et al 2000). The measured values of ro for 
APM clusters are in good agreement with the predictions of 
ACDM. They lie significantly above the rCDM predictions 
for the smallest values of d c . For the R > 1 Abell clusters, 
with d c = 52ft -1 Mpc, Peacock & West (1992) estimated 
ro = 21.1 ± 1.3, which is close to the ACDM predictions 
- 18.5/i _1 Mpc from the simulation, or 20/i~ 1 Mpc from the 
analytic theory - and also agrees with the APM results on 
this scale. 



4 POWER SPECTRA FOR THE CLUSTER 
DISTRIBUTION 

We have computed the power spectra for the cluster dis- 
tribution in our two simulations. As a comparison observa- 
tional sample we take the APM clusters analyzed by Tadros 
et al. (1998). The number density in this sample is 3.4 x 
10" 5 ih' 1 Mpc)~ 3 which is equivalent to d c = 30.9 h~ x Mpc. 
At this separation, the rCDM and ACDM simulations con- 
tain samples of about 270,000 and 915,000 clusters respec- 
tively. The upper panels of figures ^ and |s| show cluster power 
spectra from our simulations at this value of d c (filled cir- 
cles), the observational points of Tadros et al. (open trian- 
gles), and power spectra for the dark matter (open squares). 
The bias, defined as the square root of the ratio of cluster 
to dark matter power spectrum, is plotted in the lower pan- 
els. The power spectra are quite noisy at the largest scales 
because of the small number of modes in the simulated vol- 
ume. The peak in the power spectrum is nevertheless quite 
clear. The bias is nearly constant over a wide range of scales, 
and its value is close to that predicted by the SMT formulae 
(about 15% below the prediction of eqn (3)). The agreement 
is remarkable given the simplicity of the analytic arguments. 
The observed power spectra of Tadros et al. (1998) lie above 
both models by a factor of about 1.5. This is a little sur- 
prising since the correlation strength given by Croft et al. 
(1997) for the corresponding sample is quite similar to that 
predicted (see figure ||). Of course, our numerical results are 
in real space, whereas the APM power spectra are measured 
in redshift space. For these large scales, Kaiser's (1987) ex- 
pression should be applicable: 

^ = l + 2/?/3+/3 2 /5, (4) 

where (3 = Q os /b (see Eke et all996). For rCDM and 
ACDM, the correction factors are respectively 1.22 and 
1.15, less than half the observed offset between models and 
data. The remaining differences are not large and may re- 
flect residual systematics in the observational data analysis. 
Comparing the observational points with the simulations it 



appears premature to argue that a peak in the observed 
power spectrum has been detected. 



5 SUMMARY 

We have presented results for the second order clustering 
statistics of mass-limited samples of galaxy clusters in our 
Hubble Volume simulations. These simulations follow the 
matter distribution in by far the largest volumes treated to 
date, and as a result we are able to estimate clustering statis- 
tics with unprecedented precision. The two simulations we 
have studied are a rCDM universe with SI = 1 and a ACDM 
universe with Q — 0.3. Both are consistent with the fluctu- 
ation amplitude measured by COBE and with the observed 
abundance of rich clusters at z = 0. Cluster correlations are 
very similar in these two models, although slightly stronger 
in the low density case. In both cosmologies, the correlation 
length of rich clusters increases from 12 - 13 /i _1 Mpc for rel- 
atively low mass objects with mean separation 30 /i _1 Mpc 
to 22 - 27 /i~ 1 Mpc for rarer and more massive objects with 
mean separation 100 /i~ 1 Mpc. For both models, the power 
spectrum of the cluster distribution is accurately parallel 
to that of the dark matter for wavenumbers k = 0.01 - 
O.lh Mpc" 1 . 

We have compared our results with predictions from 
the analytic model of Sheth, Mo & Tormen (2000). When 
clustering strengths are compared as a function of the mean 
separation of the cluster sample, there are no adjustable pa- 
rameters and it is thus remarkable that we find good agree- 
ment in all cases. Correlation lengths are predicted by the 
analytic model to better than 10%, and the mean bias of 
the power spectrum is predicted even more accurately on 
the scales most relevant for real samples. 

We have also compared our results with published data 
on the APM cluster sample (Croft et al. 1997; Tadros et 
al. 1998). The observed trend of clustering with richness is 
very similar to those predicted in our CDM models. The 
observed correlation lengths are consistent with those pre- 
dicted by our ACDM model at all richness levels, and are 
also compatible with our tCDM model except perhaps for 
the poorest systems. The published power spectrum for the 
APM clusters agrees in shape with that predicted by the two 
models, but its amplitude is greater by about 50%. Only part 
of this discrepancy can be attributed to redshift distortion 
effects. Since the observed spectrum is based on only 364 
clusters, it may be prudent to wait for larger samples before 
drawing substantive conclusions from this disagreement. 
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